function Ef = FindFermi(x,Vx,BGCharge,E,Se1,Se2,Ef0)
% Ef0 is the initial guess for fermi level

% Grid Setup
he = E(2)-E(1);
hx = x(2)-E(1);

% Starting from Ef0, if TotalCharge>0, Ef goes up; if TotalChargae<0, Ef
% goes down
ne = 0;
Ef = Ef0+ne*he;
S = sum(BGCharge-V2Charge(Se1,Se2,E,Ef+Vx));
if S>0
    while S>0
        ne = ne+1;
        Ef = Ef0+ne*he;
        Sband = V2Charge(Se1,Se2,E,Ef+Vx);
        S = sum(BGCharge-Sband);
    end 
    S_prev = sum(BGCharge-V2Charge(Se1,Se2,E,Ef+Vx-he));
    if abs(S_prev) < abs(S)
        ne = ne-1;
    end
else
    while S<0
        ne = ne-1;
        Ef = Ef0+ne*he;
        S = sum(BGCharge-V2Charge(Se1,Se2,E,Ef+Vx));
    end
    S_prev = sum(BGCharge-V2Charge(Se1,Se2,E,Ef+Vx+he));
    if abs(S_prev) < abs(S)
        ne = ne+1;
    end
end
Ef = Ef0+ne*he;

end